In this document, I will show our main models and their respective diagnostics plots. The full explanation about these analyses can be found in our preregistered analysis plan: https://osf.io/kvuqr/

In case you are interested in the R code, you can click in the “Code” → “Show All Code” option in the right upper corner of this document.

pacman::p_load(tidyverse, # Data wrangling
               brms, # to fit Bayesian models
               tidybayes, # to wrangle and plot brms data
               rio, # to import/export files
               here, # reproducible file paths
               metafor) # to calculate log OR

# Load function to plot diagnostics

source(here("final_analyses/functions/diag_plot.R"))

set.seed(123)
# Load original data file
d = import(here("final_analyses", "data", "mortality_data.xlsx"))

# Remove studies that have 0 total events in both treatment arms (3 studies)
d = 
  d %>% 
  filter((control_events + trt_events) != 0)

# Calculate log odds ratio

d_logOR = 
  escalc(
  measure = "OR", # log odds ratio,
  
  # Tocilizumab/Sarilumab
  ai = trt_events,
  n1i = trt_total,
  
  # Control
  ci = control_events,
  n2i = control_total,
  
  data = d
) %>%
  as_tibble() %>% 
  # Dummy approach, tocilizumab = 0, sarilumab = 1
  mutate(treatment = ifelse(treatment == "tocilizumab", 0, 1)) %>% 
  # Calculate standard error
  mutate(sei = sqrt(vi)) 

Sarilumab only model

Model

\[ \begin{align*} y_i & \sim Normal(\theta_i, \sigma_i^2) \tag{Likelihood}\\ \theta_i & \sim Normal(\mu_{sarilumab}, \tau^2)\\ \\ \mu_{sarilumab} & \sim \operatorname{Normal}(0, 1.5^2) \tag{Priors}\\ \tau & \sim \operatorname{Half-Normal}(0.5) \\ \end{align*} \]

## Formula
mf = 
  # https://bookdown.org/content/3890/horoscopes-insights.html#use-the-0-intercept-syntax
  formula(yi | se(sei) ~ 0 + Intercept + (1 | study))

## Priors
priors = 
  prior(normal(0, 1.5), class = "b", coef = "Intercept") + 
  prior(normal(0, 0.5), class = "sd") # Half-Normal(0.5)

## Fit model (suppressed)

# model_sarilumab_only = 
#   brm(
#   data = d_logOR %>% filter(treatment == "1"), # only sarilumab
#   family = gaussian,
#   
#   formula = mf,
#   prior = priors,
#   sample_prior = TRUE,
#   
#   
#   backend = "cmdstanr", # faster
#   cores = parallel::detectCores(),
#   chains = 4,
#   warmup = 2000,
#   iter = 4000,
#   control = list(adapt_delta = .95),
#   seed = 123,
#   
#   file = here("final_analyses/output/fits/main/model_sarilumab_only")
# )

# Load model
model_sarilumab_only = 
  readRDS(here("final_analyses/output/fits/main/model_sarilumab_only.rds"))
print(model_sarilumab_only)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: yi | se(sei) ~ 0 + Intercept + (1 | study) 
##    Data: d_logOR %>% filter(treatment == "1") (Number of observations: 9) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~study (Number of levels: 9) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.23      0.17     0.01     0.65 1.00     2875     2635
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     0.02      0.18    -0.38     0.33 1.00     4396     3463
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.00      0.00     0.00     0.00 1.00     8000     8000
## 
## Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Diagnostics plots

diag_plot(model = model_sarilumab_only,
          pars_list = c("b_Intercept", "sd_study__Intercept"),
          ncol_trace = 4)

Posterior predictive check

pp_check(model_sarilumab_only, nsamples = 50)

Sarilumab vs. Tocilizumab model

Model

\[ \begin{align*} y_i & \sim Normal(\theta_i, \sigma_i^2) \tag{Likelihood} \\ \theta_i & \sim Normal(\mu, \tau^2)\\ \mu &= \beta_0 + \beta_1 x\\ \\ \beta_0 & \sim \operatorname{Normal}(0, 1.5^2) \tag{Priors} \\ \beta_1 & \sim \operatorname{Normal}(0, 1^2) \\ \tau & \sim \operatorname{Half-Normal}(0.5) \\ \end{align*} \]

# Formula
mf = 
  formula(yi | se(sei) ~ 0 + Intercept + treatment + (1 | study))

priors = 
  prior(normal(0, 1.5), class = "b", coef = "Intercept") +
  prior(normal(0, 1), class = "b", coef = "treatment") +
  prior(normal(0, 0.5), class = "sd") 


# model_sarilumab_vs_tocilizumab =
#   brm(
#   data = d_logOR,
#   family = gaussian,
#   
#   formula = mf,
#   prior = priors,
#   sample_prior = TRUE,
#   
#   
#   backend = "cmdstanr", # faster
#   cores = parallel::detectCores(),
#   chains = 4,
#   warmup = 2000,
#   iter = 4000,
#   control = list(adapt_delta = .95),
#   seed = 123,
#   
#   file = here("final_analyses/output/fits/main/model_sarilumab_vs_tocilizumab")
# )

model_sarilumab_vs_tocilizumab =
  readRDS(here("final_analyses/output/fits/main/model_sarilumab_vs_tocilizumab.rds"))
print(model_sarilumab_vs_tocilizumab)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: yi | se(sei) ~ 0 + Intercept + treatment + (1 | study) 
##    Data: d_logOR (Number of observations: 25) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~study (Number of levels: 24) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.13      0.09     0.01     0.35 1.00     2467     3454
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.18      0.09    -0.34     0.00 1.00     3701     3585
## treatment     0.23      0.16    -0.11     0.52 1.00     5149     4182
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.00      0.00     0.00     0.00 1.00     8000     8000
## 
## Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Diagnostics plots

diag_plot(model = model_sarilumab_vs_tocilizumab,
          pars_list = c("b_Intercept", "b_treatment", "sd_study__Intercept"),
          ncol_trace = 4)

Posterior predictive check

pp_check(model_sarilumab_vs_tocilizumab, nsamples = 50)

Sarilumab + Tocilizumab as prior model

Tocilizumab only model

\[ \begin{equation*} \sigma_{i[W]}^2 = \frac{1}{Wa_i+\frac{1}{2}} + \frac{1}{Wb_i+\frac{1}{2}} + \frac{1}{Wc_i+\frac{1}{2}} + \frac{1}{Wd_i+\frac{1}{2}} \end{equation*} \]

where \(a_i\), \(b_i\), \(c_i\) and \(d_i\) represent the number of events (death and no death) in a tocilizumab study \(i\).

W_fun = function(x) {
  sew = sqrt((1/(W*x$a + 1/2) + 1/(W*x$b + 1/2) + 1/(W*x$c + 1/2) + 1/(W*x$d + 1/2)))
  sew
}
# Create data frame with only relevant data
toci_data = d_logOR %>% 
  filter(treatment == 0) %>% 
  mutate(a = trt_events,
         b = control_events,
         c = trt_total - trt_events,
         d = control_total - control_events) %>% 
  select(study, yi, sei, a:d) %>%
  as.data.frame()

# Vector of weights, ranging from 0.01 to 1

W = c(0.01, 0.05, seq(from = 0.1, to = 1, by = 0.1))

# Apply W_fun to create a data frame where each new column is the standard error
# per W
weight_data = plyr::adply(.data = toci_data, 
            .margins = 1, 
            .fun = function(x) W_fun(x)) %>% 
  
  # Select revelant columns
  select(study, yi, V1:last_col()) %>% 
  # Change new columns names
  setNames(c(colnames(toci_data)[1:2], paste0("seW_", W)))

Now, let’s fit several frequentist random-effect meta-analysis (one per weight)

\[ \begin{align*} y_i & \sim Normal(\theta_{i[W]}, \sigma_{i[W]}^2) \tag{Likelihood}\\ \theta_{i[W]} & \sim Normal(\mu_{tocilizumab[W]}, \tau_{[W]}^2)\\ \end{align*} \]

# Create vector with all standard error column names
column_names =
  colnames(weight_data) %>% 
  data.frame() %>% 
  rename(weights = ".") %>% 
  slice(3:n())

# Create empty list to fill below with all fits in the loop
tocilizumab_freq_models_w = list()

# Create empty data frame to fill with overall mean and sd from each fit per
# weight
dat = data.frame(weight = c(0.01, 0.05, seq(from = 0.1, to =  1, by = 0.1)),
                  mean = NA, sd = NA)

# Run loop to fit all 12 fits and save info in data frame and list
for (i in 1:nrow(column_names)) {
  
  # Pick a column name
  weight_value = column_names %>% slice(i) %>% pull()
  
  # Transform data object into the long format to be able to use filter() below
  d_w = 
    weight_data %>% 
  pivot_longer(seW_0.01:last_col(),
               names_to = "sei")
  
  # Filter only data regarding the weight_value above
  d_rma = 
    d_w %>% 
  filter(sei == weight_value)
  
  # Fit meta-analysis
  fit = rma(yi = yi, sei = value,
            data = d_rma,
            method = "REML", slab = study)
  
  # Save all fits in this list
  tocilizumab_freq_models_w[[i]] = fit
  
  # Save overall mean ($beta) and sd ($se) per model in this data frame
  dat[i, "mean"] = tocilizumab_freq_models_w[[i]]$beta
  dat[i, "sd"] = tocilizumab_freq_models_w[[i]]$se
  
}

## Save output
# saveRDS(tocilizumab_freq_models_w,
#         here("final_analyses/output/fits/main/frequentist_models_prior_fits.rds"))
# 
# saveRDS(dat,
#         here("final_analyses/output/fits/main/frequentist_models_prior_mean_sd.rds"))

Sarilumab (likelihood) + Tocilizumab (prior) model

\[ \begin{align*} y_i & \sim Normal(\theta_i, \sigma_i^2) \tag{Likelihood}\\ \theta_i & \sim Normal(\mu_{sarilumab[W]}, \tau^2)\\ \\ \mu_{sarilumab{[W]}} & = \mu_{tocilizumab[W]} \tag{Priors} \\ \tau & \sim \operatorname{Half-Normal}(0.5) \end{align*} \]

# Same prior for tau in all models
tau_prior = prior_string("normal(0, 0.5)", class = "sd")

# Different priors for mu_sarilumab. 1 per W
# Define priors with string from arguments saved in dat (mean and sd)

# Weight = 0.01

## Create string to input it in prior_string() below
str_001 = paste0("normal(",
                 dat[1,2], # Mean
                 ",",
                 dat[1,3], # SD
                 ")")

## Create prior based on string above
mu_prior_001 = prior_string(str_001, class = "b", coef = "Intercept")

# Weight = 0.05
str_005 = paste0("normal(",dat[2,2],",",dat[2,3],")")
mu_prior_005 = prior_string(str_005, class = "b", coef = "Intercept")

# Weight = 0.10
str_010 = paste0("normal(",dat[3,2],",",dat[3,3],")")
mu_prior_010 = prior_string(str_010, class = "b", coef = "Intercept")

# Weight = 0.20
str_020 = paste0("normal(",dat[4,2],",",dat[4,3],")")
mu_prior_020 = prior_string(str_020, class = "b", coef = "Intercept")

# Weight = 0.30
str_030 = paste0("normal(",dat[5,2],",",dat[5,3],")")
mu_prior_030 = prior_string(str_030, class = "b", coef = "Intercept")

# Weight = 0.40
str_040 = paste0("normal(",dat[6,2],",",dat[6,3],")")
mu_prior_040 = prior_string(str_040, class = "b", coef = "Intercept")

# Weight = 0.50
str_050 = paste0("normal(",dat[7,2],",",dat[7,3],")")
mu_prior_050 = prior_string(str_050, class = "b", coef = "Intercept")

# Weight = 0.60
str_060 = paste0("normal(",dat[8,2],",",dat[8,3],")")
mu_prior_060 = prior_string(str_060, class = "b", coef = "Intercept")

# Weight = 0.70
str_070 = paste0("normal(",dat[9,2],",",dat[9,3],")")
mu_prior_070 = prior_string(str_070, class = "b", coef = "Intercept")

# Weight = 0.80
str_080 = paste0("normal(",dat[10,2],",",dat[10,3],")")
mu_prior_080 = prior_string(str_080, class = "b", coef = "Intercept")

# Weight = 0.90
str_090 = paste0("normal(",dat[11,2],",",dat[11,3],")")
mu_prior_090 = prior_string(str_090, class = "b", coef = "Intercept")

# Weight = 1
str_100 = paste0("normal(",dat[12,2],",",dat[12,3],")")
mu_prior_100 = prior_string(str_100, class = "b", coef = "Intercept")

# Put all priors together
priors = list(c(mu_prior_001, tau_prior), # W = 0.01
              c(mu_prior_005, tau_prior), # W = 0.05
              c(mu_prior_010, tau_prior), # W = 0.10
              c(mu_prior_020, tau_prior), # W = 0.20
              c(mu_prior_030, tau_prior), # W = 0.30
              c(mu_prior_040, tau_prior), # W = 0.40
              c(mu_prior_050, tau_prior), # W = 0.50
              c(mu_prior_060, tau_prior), # W = 0.60
              c(mu_prior_070, tau_prior), # W = 0.70
              c(mu_prior_080, tau_prior), # W = 0.80
              c(mu_prior_090, tau_prior), # W = 0.90
              c(mu_prior_100, tau_prior)) # W = 1.00
 
# Only data from sarilumab
sari_data = d_logOR %>% filter(treatment == 1)

# Define formula
mf = 
  formula(yi | se(sei) ~ 0 + Intercept + (1 | study))

# Define arguments for brm()
base_args = list(data = sari_data,
                 family = gaussian,
                 formula = mf,
                 iter = 4000, warmup = 2000, chains = 4,
                 cores = parallel::detectCores(),
                 control = list(adapt_delta = .99),
                 sample_prior = TRUE, seed = 123,
                 backend = "cmdstanr")

# Create vector with all weights to use later
vector_weights = data.frame(weight = c(0.01, 0.05, seq(from = 0.1, to =  1, by = 0.1)))

# Create list to store all fits
models_sarilumab_w = list()

## Run loop
# for (i in 1:nrow(vector_weights)) {
#   
#   # Show what model is running
#   print(vector_weights[i, c("weight")])
#   
#   # Put respective prior together with other arguments 
#   args = c(list(prior = priors[[i]]), base_args)
#   
#   # Fit model using brm()
#   fit = do.call(brm, args)
#   
#   # Store the whole model in a list, where
#   # [[1]] = 0.01 weight, [[2]] = 0.05, [[3]] = 0.10, [[4]] = 0.20, ...
#   # [[12]] = 1
#   models_sarilumab_w[[i]] = fit
# }

## Change names of list indexing to facilitate data wrangling later
# names(models_sarilumab_w) = vector_weights$weight

## Save list with all models
# saveRDS(models_sarilumab_w, 
#         here("final_analyses/output/fits/main/models_sarilumab_w.rds"))

# Load models
models_sarilumab_w = readRDS(
  here("final_analyses/output/fits/main/models_sarilumab_w.rds"))

Priors

vector_weights = data.frame(weight = c(0.01, 0.05, seq(from = 0.1, to =  1, by = 0.1)))

for (i in 1:nrow(column_names)) {
  
  lab = paste0("W = ", vector_weights %>% slice(i) %>% pull() )
  
  forest(tocilizumab_freq_models_w[[i]], transf = exp, xlab = "Odds Ratio")
  
  grid::grid.text(lab, .65, .7, gp=grid::gpar(cex=2))
  
}

These are the overall estimates from the frequentist models

dat %>% 
  ggplot(aes(y = "", dist = distributional::dist_normal(mean, sd),
             color = sd)) +
  stat_dist_slab(fill = NA) +
  coord_cartesian(expand = FALSE) +
  labs(x = "log odds ratio", y = "Density\n") +
  theme_minimal()

dat %>% 
  ggplot(aes(y = sd, dist = distributional::dist_normal(mean, sd))) +
  stat_dist_halfeye(position = "dodge",
                    .width = .95) +
  labs(x = "log odds ratio", y = "Standard deviation of each distribution\n") +
  theme_minimal()

Diagnostics plots

Sarilumab + Tocilizumab model
# Run loop
for (i in 1:nrow(vector_weights)) {
  # Run custom function diag_plot, 1 per weight
p = diag_plot(model = models_sarilumab_w[[i]],
          pars_list = c("b_Intercept", "sd_study__Intercept"),
          ncol_trace = 4)
# Display plot
print(p)  

# Add legend to show respective weight
lab = paste0("W = ", vector_weights %>% slice(i) %>% pull() )
grid::grid.text(lab, .75, .05, gp=grid::gpar(cex=2))
}

Posterior predictive check

# Run lopp
for (i in 1:nrow(vector_weights)) {
p = pp_check(models_sarilumab_w[[i]], nsamples = 50)

print(p)

lab = paste0("W = ", vector_weights %>% slice(i) %>% pull() )
grid::grid.text(lab, .8, .9, gp=grid::gpar(cex=2))
}

sessionInfo()
## R version 4.0.4 (2021-02-15)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] bayesplot_1.8.0 cowplot_1.1.1   metafor_3.0-2   Matrix_1.3-4   
##  [5] here_1.0.1      rio_0.5.26      tidybayes_2.3.1 brms_2.15.0    
##  [9] Rcpp_1.0.7      forcats_0.5.1   stringr_1.4.0   dplyr_1.0.7    
## [13] purrr_0.3.4     readr_1.4.0     tidyr_1.1.3     tibble_3.1.2   
## [17] ggplot2_3.3.5   tidyverse_1.3.1
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.3.1          backports_1.2.1       plyr_1.8.6           
##   [4] igraph_1.2.6          svUnit_1.0.6          splines_4.0.4        
##   [7] crosstalk_1.1.1       TH.data_1.0-10        rstantools_2.1.1     
##  [10] inline_0.3.19         digest_0.6.27         htmltools_0.5.1.1    
##  [13] rsconnect_0.8.18      fansi_0.5.0           magrittr_2.0.1       
##  [16] memoise_2.0.0         openxlsx_4.2.3        modelr_0.1.8         
##  [19] RcppParallel_5.1.4    matrixStats_0.59.0    xts_0.12.1           
##  [22] sandwich_3.0-1        prettyunits_1.1.1     colorspace_2.0-2     
##  [25] rvest_1.0.0           ggdist_2.4.0          haven_2.4.1          
##  [28] xfun_0.23             callr_3.7.0           crayon_1.4.1         
##  [31] jsonlite_1.7.2        lme4_1.1-27           survival_3.2-11      
##  [34] zoo_1.8-9             glue_1.4.2            gtable_0.3.0         
##  [37] emmeans_1.6.1         V8_3.4.2              distributional_0.2.2 
##  [40] pkgbuild_1.2.0        rstan_2.21.2          abind_1.4-5          
##  [43] scales_1.1.1          mvtnorm_1.1-1         DBI_1.1.1            
##  [46] miniUI_0.1.1.1        xtable_1.8-4          foreign_0.8-81       
##  [49] stats4_4.0.4          StanHeaders_2.21.0-7  DT_0.18              
##  [52] htmlwidgets_1.5.3     httr_1.4.2            threejs_0.3.3        
##  [55] arrayhelpers_1.1-0    ellipsis_0.3.2        farver_2.1.0         
##  [58] pkgconfig_2.0.3       loo_2.4.1             sass_0.4.0           
##  [61] dbplyr_2.1.1          utf8_1.2.1            labeling_0.4.2       
##  [64] tidyselect_1.1.1      rlang_0.4.11          reshape2_1.4.4       
##  [67] later_1.2.0           munsell_0.5.0         cellranger_1.1.0     
##  [70] tools_4.0.4           cachem_1.0.5          cli_3.0.0            
##  [73] generics_0.1.0        pacman_0.5.1          mathjaxr_1.4-0       
##  [76] broom_0.7.6           ggridges_0.5.3        evaluate_0.14        
##  [79] fastmap_1.1.0         yaml_2.2.1            processx_3.5.2       
##  [82] knitr_1.33            fs_1.5.0              zip_2.2.0            
##  [85] nlme_3.1-152          mime_0.10             projpred_2.0.2       
##  [88] xml2_1.3.2            compiler_4.0.4        shinythemes_1.2.0    
##  [91] rstudioapi_0.13       curl_4.3.2            gamm4_0.2-6          
##  [94] reprex_2.0.0          bslib_0.2.5.1         stringi_1.6.2        
##  [97] highr_0.9             ps_1.6.0              Brobdingnag_1.2-6    
## [100] lattice_0.20-44       nloptr_1.2.2.2        markdown_1.1         
## [103] shinyjs_2.0.0         vctrs_0.3.8           pillar_1.6.1         
## [106] lifecycle_1.0.0       jquerylib_0.1.4       bridgesampling_1.1-2 
## [109] estimability_1.3      data.table_1.14.0     httpuv_1.6.1         
## [112] R6_2.5.0              promises_1.2.0.1      gridExtra_2.3        
## [115] codetools_0.2-18      boot_1.3-28           colourpicker_1.1.0   
## [118] MASS_7.3-54           gtools_3.8.2          assertthat_0.2.1     
## [121] rprojroot_2.0.2       withr_2.4.2           rdocsyntax_0.4.1.9000
## [124] shinystan_2.5.0       multcomp_1.4-17       mgcv_1.8-36          
## [127] parallel_4.0.4        hms_1.1.0             grid_4.0.4           
## [130] coda_0.19-4           minqa_1.2.4           cmdstanr_0.4.0       
## [133] rmarkdown_2.8         shiny_1.6.0           lubridate_1.7.10     
## [136] base64enc_0.1-3       dygraphs_1.1.1.6